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Abstract: 


This  memorandum  presents  a  rapid  and  accurate  solution  for  the  Poisson 
equation  with  Neumann  boundary  conditions  on  general  domains.  The  stream- 
like  function  formulation  is  used  to  obtain  two  modified  Poisson  equations 
with  Dlrichlet  boundary  conditions.  The  first  equation  is  solved  for  the 
flux  components  of  the  dependent  variable,  and  are  used  in  calculating 
Dirlchlet  boundary  conditions  for  the  second  Poisson  equation.  Two  integral 
constraints  are  to  be  satisfied  in  the  present  formulation.  The  first 
constraint  is  a  consequence  of  Green's  first  integral  theorem  and  the  *- 
second  constraint  is  a  result  of  the  present  analysis.  A  new  procedure  is'^ 
implemented  in  the  boundary  conditions  to  satisfy  the  second  integral 
constraint.  This  technique  is  applied  to  calculate  the  static  pressure 
from  a  Poisson  equation  with  Neumann  boundary  conditions  in  the  solution  of 
Navler-Stokes  equations  (velocity  pressure  formulation) .  Numeridsl  results 
for  the  incompressible  viscous  flow  in  a  driven  cavity  are  presentpd  and 
compared  with  the  available  numerical  results  in  the  literature. 
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Nomenclature 

area  of  the  solution  domain 
contour,  enclosing  the  area,  A 
arbitrary  constants 
parameter  defined  In  Eq.  (17) 

Incremental  distance  along  contour,  C 
Neumann  boundary  conditions 

outward  normal  to  the  boundaries  of  the  solution  domain 
dependent  variable  and  static  pressure 
Reynolds  number 

flux  components  defined  In  Eqs.  (4)  and  (5) 
velocity  component  in  the  x-dlrectlon 
velocity  component  in  the  y-direction 
arbitrary  points  on  the  Integration  reference  line 
source  term,  Eqs.  (1)  and  (17) 


X 


streamlike  function 
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IMTHODUCTION 


A  rapid  and  accurate  solution  to  the  discrete  Poisson  equation  Is 
essential  for  the  computational  fluid  dynamic  applications  and  other 
fields  [1-8].*  In  recent  years,  highly  efficient  direct  solutions  for 
the  Poisson  equation  have  been  developed  using  cyclic  reduction  methods, 
tensor  product  methods,  the  Fourier  series  method,  and  a  few  others 
[9-13].  All  these  methods  have  one  or  more  limitations  [1]. 

The  existence  and  uniqueness  of  a  solution  to  the  Poisson  equation 
with  Neumann  (derivative)  boundary  conditions  requires  an  extra  constraint 
which  relates  the  source  term  and  the  Neumann  conditions.  This  condition 
Is  not  exactly  satisfied  In  the  numerical  solutions  due  to  truncation 
errors  which  lead  to  slow  convergence.  The  slow  convergence  of  the 
solution  for  the  Poisson  equation  has  restricted  many  of  the  numerical 
solutions  [1].  It  Is  worth  mentioning  here  that  two  methods  to  Improve 
the  convergence  of  the  numerical  techniques  are  given  In  References  [2] 
and  [3].  Briley  [2]  modified  the  source  term  In  the  Poisson  equation  to 
match  the  boundary  conditions,  while  Mlyakoda  [3]  modified  the  boundary 
conditions  Instead. 

This  memorandum  presents  a  new  method  for  solving  the  Poisson 
equation  with  Neumann  boundary  conditions  on  general  domains.  The 
streamllke  function  formulation  [15]  Is  used  to  transform  the  Poisson 
equation  with  Neumann  boundary  conditions  Into  two  Poisson  equations 
with  Dlrlchlet  boundary  conditions.  These  two  Poisson  equations  can  be 
solved  numerically  using  Iterative  and  Direct  techniques  [9-14].  The 
Direct  solvers  are  applicable  In  the  present  case  because  of  the 
resulting  Dlrlchlet  boundary  conditions. 

This  technique  Is  applied  to  calculate  the  static  pressure  from  a 
Poisson  equation  with  Neumann  boundary  conditions  In  the  solution  of 
Navler-Stokes  equations  (velocity  pressure  formulation) .  Numerical 
results  are  obtained  for  the  Incompressible  viscous  flow  In  a  driven 
cavity.  The  momentum  equations  are  Integrated  for  the  velocity  components 
using  an  explicit  marching  technique,  while  the  two  Poisson  equations  for 
the  pressure  are  solved  using  the  successive  over-relaxation  method.  The 
computed  results  are  presented  and  compared  with  the  numerical  results  of 
References  [16-20]. 

MATHEMATICAL  FORMULATION 


The  two-dimensional  Poisson  equation  with  Neumann  (derivative) 
boundary  conditions  can  be  written  In  the  following  form:** 


*Numbers  In  brackets  [  ]  designate  References,  Page  12. 

**Thls  analysis  Is  performed  for  the  two-dimensional  Poisson  equation, 
but  It  Is  directly  valid  for  the  axlsymmetrlc  case. 
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and 


0(x,y) 


(1) 


I  -  (2) 

where  F  Is  the  dependent  variable,  a  Is  the  source  term,  n  is  the  outward 
oomal  to  the  boundaries  of  the  solution  domain,  and  £  is  the  derivative 
of  the  dependent  variable  specified  along  the  boundary  contour  C. 

The  solution  to  Eq.  (1)  with  boundary  conditions,  Eq.  (2),  exists  if 
the  following  Integral  constraint  which  results  from  Green's  first  Integral 
theorem  is  satisfied: 


o(x,y)  dA 


f(C)  dC 


(3) 


where  A  is  the  area  of  the  solution  domain  and  dC  is  the  incremental  distance 
along  C. 

The  solution  to  Eq.  (1)  is  obtained  using  the  streamlike  function 
formulation  [15]  as  follows: 


Let 


and 


Q 


3y 


(A) 


(5) 


Using  the  above  equations,  Eq.  (1)  can  be  written  in  the  following 
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Equations  (6)  and  (7)  are  solved  for  the  flux  cos^neats  R  and  Q 
using  the  streamlike  function  formulation  [15].  These  ti^.,coBq>onents 
are  expressed  In  terms  of  the  streamlike  function,  x*  aod  source 
term,  o,  in  Eq.  (6)  as:  ..  J 


R 


A 

1^  +  I  o(x,y)  dx 

*1 


(8 


and 


3x 


(S, 


Equations  (8)  and  (9)  satisfy  Eq.  (6)  Identically  and  transform 
Eq.  (7)  into  the  following  second-order  differential  equation  in  x: 


IJL  +  lJL 
3x^  3y^ 


a(x,y)  dx 


(10) 


The  boundary  conditions  for  Eq.  (10)  in  terms  of  the  streamlike 
function,  x>  obtained  from  the  Neumann  conditions  for  P  given  by 
Eq.  (2). 


X 


I 


f(C)  dC  -  j  j  a(x,y)  dx  dy  + 
*1 


(11) 


where  xj^  and  Sj^  are  arbitrary  points  on  the  integration  reference  line 
and  is  an  arbitrary  constant. 


From  the  solution  of  Eq.  (10)  with  the  boundary  conditions,  Eq.  (11), 
the  streamlike  function  is  known  everywhere  in  the  solution  domain.  The 
flux  components  R  and  Q  are  calculated  from  the  streamlike  function,  x> 
using  Eqs.  (8)  and  (9).  The  integration  of  these  flux  coiq>onents  along 
the  boundaries  will  result  in  the  following  Dirlchlet  boundary  condition 
for  P: 


X  XX 

P-|  l^dx  +  l  I  a(x,y)  dx  dx  -  I  IJ  dy  +  Cj  (12) 

*1  *1  *1 


where  C.  is  an  arbitrary  constant 
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The  solution  of  Eq.  (1)  with  the  Dlrlchlet  boundary  conditions, 

Eq.  (12),  for  the  dependent  variable,  ’,  closes  the  problem. 

Integral  Constraints 

Two  Integral  constraints  have  to  be  satisfied  In  the  numerical 
solution  of  Eqs.  (1)  and  (10).  These  constraints  result  from  the 
Integration  of  Eqs.  (6)  and  (7)  over  the  area  of  the  sol  tlon  domain. 

The  first  Integral  constraint  Is  the  existence  condition  given  by  Eq.  (3). 
This  condition  exists  In  the  Dlrlchlet  boundary  conditions  when  the 
Integrals  In  Eq.  (11)  are  evaluated  over  the  total  area  and  the  closed 
contour,  C.  The  presence  of  the  Integral  constraint,  Eq.  (3),  In  the 
boundary  conditions,  Eq.  (11),  causes  a  discontinuity  In  the  streamllke 
function,  x<  point  on  the  boundary  If  the  constraint,  Eq.  (3), 

Is  not  exactly  satisfied.  This  discontinuity  Is  treated  In  a  manner 
equivalent  to  the  treatment  of  the  existence  condition  In  Reference  [2]. 
The  area  average  of  the  error  between  the  two  sides  of  Eq.  (3)  Is 
subtracted  uniformly  from  the  source  term  of  Eq.  (1). 

The  second  Integral  constraint  for  the  solution  of  Eq.  (1)  Is 


I  |§  4c  -  0  .  (i; 

A  similar  technique  to  the  treatment  of  the  first  Integral  constraint, 
Eq.  (3),  was  tried  to  satisfy  Eq.  (13),  but  without  success.  In  the 
present  study, Eq.  (13)  Is  satisfied  by  using  Neumann  conditions  for  P  on 
one  portion  of  the  contour,  C.  This  procedure  Is  simple  and  showed 
excellent  results. 

APPLICATION 


The  Incompressible  viscous  flow  in  a  driven  cavity.  Figure  1,  Is 
chosen  as  a  test  problem  In  the  present  study  because  of  Its  simple 
geometry  and  the  available  numerical  results  for  comparison.  In  the 
velocity  pressure  formulation  of  the  Navier-Stokes  equations,  the  pres¬ 
sure  is  recovered  from  a  Poisson  equation  with  Neumann  boundary  condi¬ 
tions.  In  this  equation,  conservation  of  mass  is  enforced. 

Governing  Equations 

The  two-dimensional  Navier-Stokes  equations  are  given  below  in 
nondlmenslonal  form.  The  velocity  of  the  moving  boundary  and  the 
cavity  width  are  used  as  the  characteristic  parameters. 


<  ■- W.iv-H..- 
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3u  3u  3u 

3t  3x  3y 


and  the  continuity  equation 


(14) 


(15) 


3u 

3x 


+ 


(16) 


where  u  Is  the  velocity  component  In  the  x-dlrectlon,  v  Is  the  velocity 
component  In  the  y-dlrectlon,  P  Is  the  static  pressure,  and  Re  Is  the 
Reynolds  number. 

The  following  Poisson  equation  for  the  pressure  Is  obtained  by 
adding  the  derivatives  of  Eqs.  (14)  and  (15),  with  respect  to  x  and  y, 
respectively; 


3  P  ,  3  P  f  V  /•!  7-\ 

— j  +  — ^  ■  a(x,y)  (1/a) 

y 


where 


o(x,y) 


3t 


_3_ 

3y 


3 

3u  . 

3ul 

■  3x 

u  -r-  + 

[  3x 

^  3yJ 

• 

9v 

« 

1 

3^D  .  3^D 

“37 

+  V  -5— 

3y 

Re 

2  2 
3y''J 

(17b) 


and 


D 


_  3u  .  3v 

=  -5^  37 


(17c) 
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The  continuity  equation,  Eq.  (16),  is  enforced  in  the  first  term  on 
the  right-hand  side  of  Eq.  (17). 

Equations  (14)  and  (15)  are  solved  for  the  velocity  components  u  and 
V,  respectively,  while  Eq.  (17)  is  solved  for  the  pressure,  P.  Boundary 
conditions  for  the  above  equations  in  terms  of  the  primitive  variables 
are  shown  in  Figure  2  and  Dlrlchlet  boundary  conditions  for  the  streamlike 
function,  Xi  ^nd  the  pressure,  P,  are  shown  in  Figure  3.  In  this  case, 
the  integration  reference  line  is  chosen  parallel  to  the  y-axls  and 
passing  through  the  center  of  the  lower  wall.  Figure  3.  The  Neumann 
conditions  for  the  pressure  are  retained  on  the  upper  moving  wall  for 
two  reasons;  the  first  is  to  satisfy  the  Integral  constraint,  Eq.  (13)  : 
and  the  second  reason  is  to  eliminate  the  effect  of  the  singularities  at 
the  upper  corners  on  the  solution  of  the  second  Poisson  equation. 

RESULTS 


A  numerical  code  is  written  in  FORTRAN  IV  for  the  solution  of  the 
governing  Eqs.  (14),  (15),  and  (17).  The  momentum  Eqs.  (14)  and  (15) 
are  integrated  for  the  velocity  components  u  and  v,  respectively,  using 
an  explicit  marching  technique  [1].  The  static  pressure  is  obtained 
from  the  solution  of  Eq.  (17)  using  the  present  method.  The  numerical 
solutions  of  the  Poisson  equations  are  calculated  using  the  Successive 
Over-Relaxation  method  (SOR) . 

The  numerical  code  for  the  Poisson  equation  has  been  tested  against 
analytical  solutions.  Secondly,  the  numerical  solutions  of  a  fully- 
developed  flow  in  horizontal  and  vertical  channels  have  been  used  to 
test  the  combination  of  the  marching  explicit  solution  and  Poisson 
solver . 

Figures  4  and  5  show  the  stream  function  and  the  static  pressure 
contours  for  R^  «  1.  The  streamlines  in  Figure  4  are  symmetric  about 
the  cavity  vertical  centerline,  as  expected,  and  are  plotted  here  for 
a  checkup.  The  static  pressure  contours  shown  in  Figure  5  are  almost 
antisymmetric  about  the  cavity  vertical  centerline.  In  this  case,  the 
static  pressure  is  set  equal  to  zero  at  the  center  of  the  lower  wall. 

Two  groups  of  results  are  found  in  the  literature  for  the  cavity 

problem  and  are  used  for  comparison  in  the  present  study.  The  first 
set  of  results  is  obtained  using  finite  element  techniques  and  the 
second  set  of  results  is  obtained  using  finite  difference  methods. 

Figure  6  shows  the  velocity  component  u  for  R^  -  10  on  the  horizontal 

line  y  ■  1/14  from  the  moving  wall.  It  can  be  seen  from  Figure  6  that 

the  present  results  are  in  good  agreement  with  the  results  of  the  finite 
elements  [19-20]  and  the  finite  difference  methods  [16-18].  Figures  7 
8,  and  9  show  the  velocity  components  for  R^  ■  1.  The  velocity  component  u 
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on  the  vertical  line  x  *  0.1  is  shown  in  Figure  7.  The  present  results 
In  Figure  7  are  In  good  agreement  with  the  finite  difference  solution 
of  Reference  [16]  using  the  vortlclty  stream  function  formulation.  In 
Figure  7,  It  can  be  seen  that  the  primitive  variable  solutions  of 
Reference  [19]  show  wiggles  at  points  y  ■  0.7  and  y  “  0.8.  The  velocity 
component  v  for  the  same  Reynolds  number  Rg  -  1  Is  shown  In  Figure  8  on 
the  horizontal  centerline.  The  finite  element  methods  predict  higher 
values  for  the  velocity  component  v  than  the  finite  difference  solutions. 
The  same  trend  Is  shown  In  Figure  9  for  the  velocity  component  u  on  the 
vertical  centerline. 

CONCLUSIONS 


A  method  to  obtain  converged  solutions  for  the  Poisson  equation 
with  Neumann  boundary  conditions  Is  presented.  The  computation  time 
In  the  present  method  Is  approximately  twice  the  computation  time 
required  for  the  solution  of  a  Poisson  equation  with  Dlrlchlet  boundary 
conditions  on  the  same  domain.  The  present  method  extends  the  validity 
of  the  Direct  Solvers  to  the  solution  of  the  Poisson  equation  with 
Neumann  boundary  conditions  on  general  domains.  The  application  of 
this  method  In  the  solution  of  Navler-Stokes  equations  makes  the 
primitive  variables  formulation  competitive  with  the  vortlclty  stream 
function  formulation. 
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Figure  2.  Boundary  Conditions  In  Prialtlve  Variables. 
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Figure  3.  Dirichlet  Boundary  Conditions  for  x 
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Figure  6.  Velocity  Component  u  at  y  >1/14  from  moving  wall.  Re  ■  10 


Figure  7.  Velocity  Component  uatx-0.1,  Re-1 
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Figure  8,  Velocity  Component  v  at  y  ■  1/2,  Re  ■  1. 
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